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Abstract 

We present an analytical description of the distribution of diagonal lines in Recur- 
rence Plots (RPs) for white noise and chaotic systems, and find that the latter one 
is linked to the correlation entropy. Further we identify two scaling regions in the 
distribution of diagonals for oscillatory chaotic systems that are hinged to two pre- 
diction horizons and to the geometry of the attractor. These scaling regions cannot 
be observed with the Grassberger-Procaccia algorithm. Finally, we propose methods 
to estimate dynamical invariants from RPs. 



1 Introduction 



The Recurrence Quantification Analysis (RQA) quantifies structures found 
in Recurrence Plots (RPs) to yield a deeper understanding of the underly- 
ing process of a given time series [3,5]. Even though this method is widely 
applied [14,15,26,27,28,11,13], the scarce mathematical description is a main 
drawback. First steps in the direction of an analytical description were made 
by Faure et al. [4], who gave analytical results for the cumulative distribution 
of diagonals -P e c (/) in RPs in the case of chaotic maps and linked the slope 
of this distribution to the Kolmogorov-Sinai entropy. Gao and Cai [5] related 
the distribution -P e c (Z) to the largest Lyapunov exponent and the information 
dimension. 

In this paper we give an analytical expression for the distribution of diago- 
nals in RP in the case of stochastic processes and extend the results of [4,5] 
to chaotic flows. Further we compare our approach with the Grassberger- 
Procaccia (G-P) algorithm [6] and show some advantages of the RP method 
estimating some invariants of the dynamics, such as the correlation entropy. 
One of the most remarkable differences between our approach and the G-P al- 
gorithm, is that we find two different scaling regions for chaotic flows, such as 
the Rossler system, instead of the single one obtained with the G-P algorithm. 
This new scaling region can be linked to the geometry of the attractor and 
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defines another characteristic time scale of the system. Beyond we propose 
optimized measures for the identification of relevant structures in the RP. 
The outline of this paper is as follows. In Sec. 2 we briefly introduce RPs. After 
considering in Sec. 3 the RPs of white noise, we proceed to general chaotic sys- 
tem (Sec. 4). Then, we exemplify our theoretical results for the Rossler system 
(Sec. 5) and present the two different scaling regions that characterize the sys- 
tem. Finally, we propose to estimate main characteristics of nonlinear systems 
from the RP which extends the importance of the RQA (Sec. 6). 



2 Recurrence Plots and Recurrence Quantification Analysis 

RPs were introduced to simply visualize the behavior of trajectories in phase 
space [3]. Suppose we have a dynamical system represented by the trajectory 
{Si} for i — 1, . . . , N in a <i-dimensional phase space. Then we compute the 
matrix 

Ri tj = Q(e-\xi-Xj\), i,j = l...N, (1) 

where e is a predefined threshold and G(-) is the Heaviside function 1 . The 
graphical representation of Rjj called Recurrence Plot is yielded encoding 
the value one as "black" and zero as "white" point. A homogeneous plot with 
mainly single points may indicate a mainly stochastic system. Paling away 
from the main diagonal may indicate a drift i.e. non-stationarity of the time 
series. A main advantage of this method is that it allows to apply it to non- 
stationary data [14]. 

To quantify the structures that are found in RPs, the Recurrence Quantifi- 
cation Analysis (RQA) was proposed [26]. There are different measures that 
can be considered in the RQA. One crucial point for these measures is the 
distribution of the lengths of the diagonal lines P e (l) that are found in the 
plot. In the case of deterministic systems the diagonal lines mean that trajec- 
tories in the phase space are close to each other on time scales that correspond 
to the lengths of the diagonals. In the next sections we show that there is a 
relationship between P e (l) and the correlation entropy. On the other hand we 
compute the distribution of diagonals for random processes to see that even 
in this case, there are some diagonals which can lead to pitfalls in the inter- 
pretation of the RQA because noise is inevitable in experimental systems. A 
more detailed discussion of this problem is given in [24]. 



The norm used in Eq. 1 is in principle arbitrary. For theoretical reasons, that we 
will present later, it is preferable to use the maximum norm. However the numeri- 
cal simulations of this paper are based on the Euclidian norm to make the results 
comparable with the literature. The theoretical results of this paper hold for both 
choices of the norm. 
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3 Results for white noise 



In this section we compute analytically the probability to find a black or 
recurrence point and the distribution of diagonals of length / in the RP in the 
case of independent noise. The probability to find a recurrence point in the 
RP is given by 

1 N 

P ^) = J^£ R - ( 2 ) 

and the probability to find a diagonal of at least length / in the RP is defined 

as 

Y n l-l 

P e^ =i^LW £ II Ki + m, j+ m, (3) 
i,j=l m=0 

where c stands for cumulative. Note that = P e (X)- 

We consider a random variable X with probability density p(x). Suppose that 
{xi} for i — 1, . . . , N is a realization of X and we are interested in the dis- 
tribution of the distances of each point to all other points of the time series. 
This can be done by computing the convolution of the density p(-) 

R(x) = p(x) * p(x). (4) 

Pb{e) is then gained by integrating R(x) over [—£,£] 

P b {e) = j £ R(x)dx = 2^ R(x)dx. (5) 

Note that P&(e) is invariant against shuffling of the data. For [0, 1] uniformly 
distributed noise, R(x) is given by 

. 1 1 — \x\ if \x\ < 1 
R{X) ^ (o elsl (6) 

and hence the probability Pb{e) for RPs and CRPs is given by 

P b (e) =2e -e 2 + &(e -1) [l - 2e + e 2 ] (7) 

For Gaussian white noise one finds Pfc(e) = erf (^:), where a is the standard 
deviation. Now it is straightforward to compute P^{1) in the in CRPs (in RPs 
only asymptotically). As the noise is independent, we obtain 

P c £ (l) = P b (e) 1 . (8) 

The probability to find a recurrence point Pb{s) is in both RPs and CRPs 
independent of the preceeding point on the diagonal (except in the main di- 
agonal). Eq. (8) shows that the probability to find a line of length / decreases 
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exponentially with /. For our example of uniformly distributed noise we get 

P £ %1) = (2e - e 2 ) 1 . (9) 
Note that in this case the exponential decay depends on e. 



4 Results for chaotic systems 



We present in this section an approach for chaotic systems. It is an extension 
of the results presented in [4] for chaotic maps and also covers general chaotic 
flows. To estimate the distribution of the diagonals in the RP, we start with 
the correlation integral [7] 

C(e) = ^lim j-j^ x {number of pairs with \xi — Xj\ < e}. (10) 

Note that the definition of -P&( £ ) coincides with the definition of the correlation 
integral 

1 N i N 

i=l i,j = 1 

This fact allows to link the known results about the correlation integral and 
the structures in RPs. 

We consider a trajectory x(t) in the basin of attraction of an attractor in the 
(i-dimensional phase space and the state of the system is measured at time 
intervals r. Let {1,2,..., M(e)} be a partition of the attractor in boxes of size 
e. Then p(ii, ii) denotes the joint probability that x(t = r) is in the box i±, 
x(t = 2r) is in the box i 2 , and x(t = It) is in the box i\. The order-2 Renyi 
entropy [19,8] is then defined as 

K 2 = — lim lim lim — In V] p 2 (h, ■ ■ ■ , k)- (12) 

i\,...,H 

We can approximate p(ii, . . . ,ii) by the probability P t j(x,e) of finding a se- 
quence of points in boxes of length e about x(t = r), x(t = 2r), x(t = It). 
Assuming that the system is ergodic, which is always the case for chaotic 
systems as they are mixing, we obtain 

■y N y N 

]T p*(ii,...,ii) = -^E^i'---'*') ~ jyE p t-i( f ' £ )> (I 3 ) 

h,...,il t=l t=l 

where p t (ii, . . . represents the probability of being in the box i\ at time 
t = t, in the box i 2 at time t = 2r, ... and in the box i\ at time t = It. Further 
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we can express Ptj(x,e) by means of the recurrence matrix 

■y N l-l ^ N l-l 

P tA%, £ ) = T7 S II ( £ ~ Wt+m - Xs+m\) = T7 H II ^t+m,s+m- (14) 
iV s=lm=0 iV s=lm=0 

Hence we obtain an estimator for the order-2 Renyi entropy by means of the 
RP 

K 2 (e,l) = lln (-L £ n Rt + m,s +m ) ■ (15) 



(*) 

Note that (*) is the cumulative distribution of diagonal lines -P £ c (/) (Eq. (3)). 
Therefore, if we represent -P e c (/) in a logarithmic scale versus / we should obtain 
a straight line with slope — K 2 {s)t for large Vs. 

On the other hand, in the G-P algorithm the /-dimensional correlation integral 
is defined as 

1 N ( f 1 ' 1 2 \ 1/2 \ 

Grassberger and Procaccia [10] state that due to the exponential divergence 
of the trajectories, requiring 

i-i 

Wi+k - Xj+kf < e 2 (17) 

A:=0 

is essentially equivalent to 

\x i+k - Xj+k\ < ^ for k = 1, . . . , I (18) 
which leads to the ansatz: 

C,(e) ~£ I/ exp(-ZrK 2 ). (19) 

Further they make use of Takens embedding theorem [23] and reconstruct the 
whole trajectory from / measurements of any single coordinate. Hence they 
consider 

Q(e) = Jjm ± £ (e - (j£ \x i+k - x j+k \ 2 ^ j (20) 

and use the same ansatz Eq. (19) for Ci(e). Then, the G-P algorithm obtains 
an estimator of K 2 considering 

K^) = \\*-2&- (21) 



■5 



Due to the similarity of the RP approach to the G-P one, we state 

P £ C (Z) ~ £ p 2 (i u ...,U) ~ Q(e) ~ e^exp(-/r^ 2 ). (22) 

The difference between both approaches is that in -P £ c (/) we further consider 
information about I vectors, whereas in Ci(e) we have just information about I 
coordinates. Besides this, in the RP approach / is a length in the plot, whereas 
in the G-P algorithm it means the embedding dimension. As K 2 is defined for 
I — > 00, the RP approach seems to be more appropriate than the G-P one, as 
it is always problematic to use very high embedding dimensions [21]. 
A further advantage of the RP method is that it does not make use of the 
approximation that Eq. (17) is essentially equivalent to Eq. (18). The quantity 
that enters the RPs is directly linked to the conditions Eq. (18) and hence uses 
one approximation less than the G-P method. 

One open question for both methods is the determination of the scaling re- 
gions. It is somewhat subjective and makes a rigorous error estimation prob- 
lematic. For the cases considered in this paper we have found that 10,000 
data points assure reliable results for both methods. Even 5,000 data points 
allow for a reasonable estimation, whereas 3,000 data points or less yield very 
small scaling regions that are difficult to identify. However, the RP method 
is advantageous for the estimation of K 2 as the representation is more direct. 
The most important advantage is presented in the next section: RPs allow to 
detect a new scaling region in the Rossler attractor that cannot be observed 
with the G-P algorithm. 



5 The Rossler System 

We analyze the Rossler system with standard parameters a = b = 0.2, c = 5.7 
[20]. We generate 15,000 data points based on the Runge Kutta method of 
fourth order and neglect the first 5,000. The integration step is h — 0.01 and 
the sampling rate is 20. 

First, we estimate K 2 by means of the G-P algorithm. Fig. 1 shows the results 
for the correlation integral in dependence on e. There is one well-expressed 
scaling region for each embedding dimension /. Then we get from the vertical 
distances between the lines an estimate of K 2 (Fig. 2), K 2 = 0.070 ± 0.003. 
Next, we calculate the cumulative distribution of the diagonal lines of the RP 
in dependence on the length of the lines / (Fig. 3). For large / and small e 
the scaling breaks down as there are not enough lines in the RP. The most 
remarkable fact in this figure is the existence of two well differentiated scaling 
regions. The first one is found for 1 < / < 84 and the second one for I > 85. 
The existence of two scaling regions is a new and striking point obtained from 
this analysis and is not observed with the G-P method. The estimate of K 2 
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Fig. 1. G-P algorithm for the Rossler system. I varies from 3 (top) to 27 (bottom) 
in steps of 3. 
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Fig. 2. Estimation of K 2 i for the Rossler system with the G-P algorithm. The line 
is plotted to guide the eye. 

from the slope of the first part of the lines is K 2 ~ 0.225 ± 0.03 (Fig. 4) and 
the one from the second part is K 2 ~ 0.0675 ± 0.004 (Fig. 5). Hence, K 2 is 
between 3-4 times higher than K 2 . As K 2 is defined for / — > oo, the second 
slope yields the estimation of the entropy. 

However, the first part of the curve is interesting too, as it is also independent 
of e. The region 1 < / < 84 characterizes the short term dynamics of the system 
up to three cycles around the fix point and corresponds in absolute units to 
a time of t = 16.8, as we use a sampling rate of St = 0.2. These three cycles 
reflect a characteristic period of the system that we will call recurrence period 
T rec . It is different from the dominant "phase period" T ph , which is given by 
the dominant frequency of the power density spectrum. T rec however, is given 
by recurrences to the same state in phase space. Recurrences are represented 
in the plot by vertical (or horizontal, as the plot is symmetric) white lines. 




7 



0° 



100 200 300 400 



Fig. 3. RP method for the Rossler system, e varies logarithmically from 10 2 to 10.0 
(bottom to top) 
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Fig. 4. RP method for the Rossler system: slope of the curves N£(l) in the first 
region for three different choices of the scaling region in I. 

Such a white line occurs at the coordinates i,j if 



R 



i,j+m 



= < 



1 if m = -l 

for me{0,...,l-l} 

1 if m = I. 



(23) 



The trajectory x n for times n = j — 1, . . . , j+l is compared to the point xV Then 
the structure given by Eq. 23 can be interpreted as follows. At time n — j — 1 
the trajectory falls within an £-box of xi. Then for n = j, . . . , j + 1 — 1 it moves 
outside of the box, until at n — j + I it recurs to the £-box of Hence, the 
length of the white line is proportional to the time that the trajectory needs 
to recur close to Xj. 

In Fig. 6 we represent the distribution of white vertical lines in the RP. The 
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Fig. 5. RP method for the Rossler system: slope of the curves N£(l) in the second 
region for three different choices of the scaling region in I 
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Fig. 6. Number of vertical white lines in the Recurrence Plot of the Rossler system 
with standard parameters, e = 0.05 and based on 60,000 data points. 

period of about 28 points corresponds to T ph . However, the highest peak is 
found at a lag of about 87 points (the second scaling region begins at I = 
85). This means that after this time most of the points recur close to their 
initial state. This time also defines the recurrence period T rec . For the Rossler 
attractor with standard parameters we find T rec = 3T ph . 
For predictions on time scales below the recurrence period, = is a 

better estimate of the prediction horizon than r = I/K2. This interesting 
result means that the possibility to predict the next value within an £-range is 
in the first part by a factor of more than 3 times worse than it is in the second 
part, i.e. there exist two time scales that characterize the attractor. The first 
slope is greater than the second one because it is more difficult to predict the 
next step if we have only information about a piece the trajectory for less 
than the recurrence period. Once we have scanned the trajectory for more 
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than T rec , the predictability increases and the slope of -P £ c (/) in the logarithmic 
plot decreases. Hence the first slope, as well as the time scale at which the 
second slope begins, reveal important characteristics of the attractor. 
To investigate how the length of the first scaling region depends on the form 
of the attractor, we have varied the parameter c of the Rossler system with 
fixed a = b = 0.1, so that different types of attractors appear [1]. Especially we 
have studied the cases c = 9, which yields T rec = 2T ph , and c = 30, which gives 
T rec = 4T ph . In both cases the length of the first scaling region corresponds as 
expected to T rec . 

On the other hand, the existence of the two scalings may be linked to the 
nonhyperbolic nature of the Rossler system for this attractor type, because 
the resulting two time scales have been also recently found by Anishchenko et 
al. based on a rather subtle method [2]. This effect also is detectable in other 
oscillating nonhyperbolic systems like the Lorenz system and will be studied 
in more detail in a forthcoming paper. 



6 Dynamical invariants for the RQA 

With regard to our theoretical findings in Sec. 4 we have to assess the quality 
of the possible results of the RQA. 

The measures considered in the RQA [26] are not invariants of the dynam- 
ical system, i.e. they usually change under coordinate transformations, and 
especially, they are in general modified by embedding [25]. Hence, we propose 
new measures to quantify the structures in the RP, that are invariants of the 
dynamical system. 

The first measure we propose, is the slope of the cumulative distribution 
of the diagonals for large /. We have seen that it is (after dividing by r) an 
estimator of the Renyi entropy of second order K 2 , which is a known invariant 
of the dynamics [22]. On the other hand, we also can consider the slope of 
the distribution for small /'s, as this slope shows a clear scaling region, too. 
The inverse of these two quantities, is then related to the forecasting time at 
different horizons. Especially the transition point from the first to the second 
scaling region is an interesting characteristic of the system. 
The second measure we introduce, is the vertical distance between P^{1) 
for different e l s. From Eq. (22) one can derive 

This is an estimator of the correlation dimension D 2 [8]. The result for the 
Rossler system is represented in Fig. 7. The mean value of D 2 (e) is in this case 
1.86 ±0.04. This result is in accordance with the estimation of D 2 by the G-P 
algorithm given in [18], where the value 1.81 is obtained. With a modified G-P 
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Fig. 7. Estimation of the correlation dimension -D2 for the Rossler attractor by the 
RP method. The parameters used for the Rossler system and the integration step 
are the same as in Sec. 5. 

algorithm a value of 1.89 was reported [18]. 

The third measure we suggest, is an estimator of the generalized mutual 
information of order 2, 

J 2 (r) = 2H 2 - H 2 (t) (25) 

where 

H 2 = -lnY.pl H 2 {r) = (26) 

i i 

are the generalized Renyi's second order entropy (also correlation entropy) 
and its corresponding joint second order entropy [17]. This measure can be 
estimated using the G-P algorithm as follows [12] 



7 2 (£,r) = ln(C 2 (£,r))-21n(C 1 ( £ )). 



(27) 



Instead, we can estimate hij) using the recurrence matrix. As discussed in 
the preceding sections, one can estimate H 2 as 



Ho = - In 



1 N 

iV ij=l 



(28) 



Analogously we can estimate the joint second order entropy by means of the 
recurrence matrix 



H 2 (t) = - In 



1 N 



(29) 



We compare the estimation of I 2 {r) based on the G-P algorithm with the one 
obtained by the RP method in Fig. 8. We see, that the RP method yields 
systematically higher estimates of the mutual information, as in the case of 
the estimation of the correlation entropy. However, the structure of the curves 
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Fig. 8. Comparison of the estimators of the mutual information for the x-component 
of the Rossler system computed by the RP method (solid line) and the G-P algorithm 
(dashed line). The parameters used for the Rossler system and the integration step 
are the same as in Sec. 5. 

is qualitatively the same (it is just shifted to higher values by about 0.2). A 
more exhaustive inspection shows, that the difference is due to the use of the 
Euclidean norm. The estimate based on the RP method is almost independent 
of the norm, whereas the estimate based on the G-P algorithm clearly depends 
on the special choice. If the maximum norm is used (in G-P and RP) both 
curves coincide. 

Note that the estimators for the invariants we propose are different from the 
ones of the G-P algorithm. Therefore, the obtained values are slightly different, 
too. 

The three measures that we have proposed, are not only applicable for chaotic 
systems but also for stochastic ones as the invariants are equally defined for 
both kinds of systems. 



7 Conclusions 

In this paper we have presented an analytical expression for the the distri- 
bution of diagonals P £ (l) for stochastic systems and chaotic flows, extending 
the results presented in [4]. We have shown that P £ (l) is linked to the 2-order 
Renyi entropy rather than to the Lyapunov exponent. Further we have found 
in the logarithmic plot of P e {l) two different scaling regions with respect to e, 
that characterize the dynamical system and are also related to the geometry 
of the attractor. This is a new point that cannot be seen by the G-P algorithm 
and will be studied in more detail in a forthcoming paper. The first scaling 
region defines a new time horizon for the description of the system for short 
time scales. Beyond the RP method does not make use of high embedding 
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dimensions, and the computational effort compared with the G-P algorithm is 
decreased. Therefore the RP method is rather advantageous than the G-P one 
for the analysis of rather small and/or noisy data sets. Besides this, we have 
proposed different measures for the RQA, like estimators of the second order 
Renyi entropy K2, the correlation dimension D2 and the mutual information, 
that are, in contrast to the usual ones, invariants of the dynamics [25]. 
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